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ABSTRACT 

Wc present a purely theoretical study of the morphological evolution of self-gravitating 
systems formed through the dissipation-less collapse of N point sources. We explore 
the effects of resolution in mass and length on the growth of triaxial structures formed 
by an instability triggered by an excess of radial orbits. We point out that as resolution 
increases, the equilibria shift, from mildly prolate, to oblate. A number of particles 
N ~ 100, 000 or larger is required for convergence of axial aspect ratios. An upper 
bound for the softening, e w 1/256, is also identified. We then study the properties 
of a set of equilibria formed from scale- free cold initial mass distributions, p oc r~'^\ 
^ 7 ^ 2. Oblateness is enhanced for initially more peaked structures (larger 7's). We 
map the run of density in space and find no evidence for a power-law inner structure 
when 7^3/2 down to a mass fraction < 0.1% of the total. However when 3/2 < 7 ^ 2 
the mass profile in equilibrium is well matched by a power-law of index « 7 out to 
a mass fraction « 10%. We interpret this in terms of less effective violent relaxation 
for more peaked profiles when more phase mixing takes place at the centre. We map 
out the velocity field of the equilibria and note that at small radii the velocity coarse- 
grained distribution function is Maxwellian to a very good approximation. 

We extend our study to non-scale-free initial conditions and finite but sub-virial 
kinetic energy. For cold collapses the equilibria are again oblate, as the scale-free 
models. With increasing kinetic energy the equilibria first shift to prolate morphology 
and then to spherical symmetry. 

Key words: numerical method: iV-body; galaxies, gravitational dynamics 



1 Introduction 

This is a purely theoretical investigation into the formation 
of equilibria reached through a phase of gravitational col- 
lapse (Lynden-Bell's [1967] violent relaxation process) using 
numerical A''-body calculations. Our motivation for this un- 
dertaking is drawn both from observations and theory, as 
discussed first below. A survey of recent work in that area 
follows. 



1.1 Cusps in triaxial galaxies 

HST images of elliptical galaxies have revealed cuspy (p oc 
r~^) density profiles down to resolution limits (e.g., Lauer 
et al. 1995; Gebhardt et al. 1996; Laine et al. 2003). Fits to 
their luminosity profiles indicate power-law indices ranging 
from 7«l/2to7«2. These observations have triggered 
several studies of the origin and orbital content of singular- 
ities at the heart of galaxies. The orbital structure of cuspy 



triaxial galaxies would harbour a broad range of resonant 
and near-resonant orbits, including also chaotic orbits (e.g., 
Merritt & Fridmann 1996; HoUey-Bockelmann et al. 2001). 
Whether or not the cusp hosts a massive (single or binary) 
black hole affects both the orbital families that support the 
cusp as well as the overall mass profile (HoUey-Bockelmann 
et al. 2002; Poon & Merritt 2001; Nakano & Makino 1999; 
Merritt & Cruz 2001). Furthermore, the diffusion of chaotic 
orbits in cuspy potentials may lead to a rapid readjustment 
of the equilibrium (Merritt & Fridman 1996; Kandrup & 
Siopis 2003). Thus the possibility of tracking individual or- 
bits in self-consistent potentials can hand diagnostics on the 
stability of ellipticals and on the demographics of black holes 
in galaxies. 

In this contribution, we construct cuspy triaxial galaxies 
obtained from the classic initial-conditions problem of vio- 
lent relaxation. In particular we are interested in the numer- 
ical resolution of central density peaks and, at large radii, 
the global morphology in equilibrium. Analytic and semi- 
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analytic models of triajcial galaxies can be constructed from 
distribution functions subject to appropriate constraints 
(see e.g., Dehnen & Gerhard 1994; HoUey-Bockelmann et al. 
2001; van de Ven et al. 2003). However, we chose the flexi- 
ble three-dimensional Af-body approach to explore a range 
of initial conditions and study the time-evolution of the sys- 
tems. This allows total freedom with regard to the symmetry 
and boundary conditions imposed on the system. Drawbacks 
include finite spatial and mass resolution imposed by limited 
computer resources. This calls on us to check which simula- 
tions parameters allow sufficient accuracy and convergence 
in the numerics, with a view to develop a full library of 
non-analytic models in the most cost-efficient way. 

We divide our study in two parts. This contribution 
is the first part devoted largely to checks of the numerical 
setup and a parameter survey. In a forthcoming contribu- 
tion we will explore the families of orbits and derive observ- 
ables from the equilibria obtained. Here, first we define the 
initial conditions for collapse in terms of an accretion prob- 
lem and give details of the numerical method (§2). Draw- 
ing from known correlations between equilibrium and initial 
energy distributions (van Albada 1982; Aguilar & Merritt 
1990; Henriksen & Widrow 1999) we setup sub-virial scale- 
free power-law mass profiles (oc r '; see Eq.[l] below). We 
monitor the convergence of physical parameters (axis ratios, 
density profiles) with particle number, A'^, and linear reso- 
lution e in §3. We then relate equilibrium properties to the 
initial conditions in terms of meiss profile and virial ratio 
(see §3 and §4), with the expectations that the power-index 
of the inner density profile will match 7 of the initial con- 
ditions. The properties of a set of cuspy density profiles are 
worked out for 7 in the range ^ 7 ^ 2, covering the range 
of power indices derived for observed ellipticals (§5 and §6). 
We conclude with an extended discussion and possible ap- 
plications of these results in a cosmological context. A brief 
survey of the literature on the topic of violent relaxation is 
helpful to set our goals in context. 

1.2 Previous numerical work on violent relaxation 

In a ground-breaking A'^-body numerical investigation of vi- 
olent relaxation, van Albada (1982) showed that the de Vau- 
couleurs Ti^^* projected luminosity profile of ellipticals could 
be understood as the outcome of gravitational collapse. This 
feature proved attractive since massive ellipticals are largely 
pressure-supported with little net rotation (Davies et al. 
1983; Binney & Merrifield 1999), a by-product of violent 
relaxation (e.g. May & van Albada 1984; Curir & Diaferio 
1994). Later, radial orbit instabilities (ROI) were shown to 
develop in numerical renditions of anisotropic systems con- 
structed from equilibrium distribution functions f{E, J^) of 
energy and square angular momentum (Merritt & Aguilar 
1985; Barnes, Hut & Goodman 1986) as well as in the end- 
results of violently relaxed systems (Aguilar & Merritt 1990; 
Canizzo & HoUister 1992). The ROI develops in e.g. collaps- 
ing (sub-virial) spherical distributions due to large radial 
bulk motion, so that the initial symmetry is lost and the 
systems become triaxial in equilibrium. The aspect ratios 
(of mean ~ 1 : 2 [E5]) attained in these studies cover the en- 
tire spectrum of ellipticals. Aguilar & Meritt (1990) demon- 
strate that the equilibrium distribution functions / satis- 
fies Antonov's stability to radial perturbations df/dE < 0, 



where E is the binding energy. However, van Albada's low- 
resolution core-halo equilibria makes direct application of 
these results to the inner structure of cuspy de Vaucouleur 
and HST galaxies less than straightforward. It was noted 
that relaxation from less smooth (clumpy) initial conditions 
leads to more compact virialised equilibria (e.g. McGlynn 
1984; Aguilar & Merritt 1990; Roy & Perez 2004) possi- 
bly offering a way forwarcQ. In the late 1980's the emphasis 
on galaxy formation shifted to include gas cooling and the 
transformation of spirals to ellipticals through mergers (see 
e.g. Barnes & Hernquist 1992). Dissipation-less dynamics 
with A''-body calculations recast in the framework of hierar- 
chical structure formation in a CDM cosmogony was found 
to give rise to equilibria with steep central cusps (7^1: 
Navarro et al. 1997, 2004; Fukushige & Makino 1997, 2001; 
Fukushige, Kawai & Makino 2004; Moore et al. 1998, 2004; 
Diemand et al. 2005). Thus repeated episodes of mass accre- 
tion would seemingly lead to equilibria with steeper central 
cusps. 

Both Aguilar & Merritt (1990) and Canizzo & HoUis- 
ter (1992) modelled sub-virial mass in-fall with power- law 
initial mass profiles distributed spherically. Each mass shell 
collapses to the origin at a different time, and the accretion 
is continued until the last mass shell converges to the origin. 
They have shown that (1) the ensuing equilibria are highly 
prolate; (2) the aspect ratio (defined from the eigenvectors 
of the inertia tensor) increases radially and is near unity at 
large radii; (3) the spherically-averaged density profiles of 
equilibrium systems correlates with the initial profile. These 
investigations were carried out using multi-polar series ex- 
pansion and the TREE (Barnes & Hut 1986) integrators re- 
spectively with N — 5, 000 and 10,000 mass elements. More 
recent work by Roy & Perez (2004) used up to A'' = 30, 000 
particles. One important result drawn from this study is 
that homogeneous initial mass profile lead to equilibria that 
match better cored globular cluster in equilibrium, a con- 
clusion that lends support to the analysis of LMC clusters 
of Boily et al. (1999). 

The current, purely theoretical study must be cast in 
the context of the interesting, full-fledged parameter survey 
of Roy & Perez to highlight differences : 

1) Recently we have shown that collapse factors in 
spherical symmetry are indistinguishable from those ob- 
tained from non-spherical in-fall when the mass resolution 
N < 10^ (Boily et al. 2002 = BAK+02). This raises is- 
sues with the growth of velocity anisotropies near the time 
of maximum contraction in studies with lower resolution. 
Aarseth et al. (1988) and Hozumi et al. (1996) have argued 
that the tangential velocity dispersion grows faster than the 
radial component at the bounce. Tangential velocities will 
arise from the growth of fragmentation modes (McGlynn 
1984; Aarseth et al. 1988). Clearly the velocity field must 
develop fully for numerical convergence of the end-product 
properties. The results of BAK-l-02 would suggest a mini- 
mum ~ 10^ in order to discriminate between spherical 
and non-spherical growth modes, by resolving the dynamics 
at maximum contraction well. We confirm in §4 that nu- 



^ The same holds if a smooth collapsing system has random tan- 
gential velocities initially (Hozumi ot al. 2000; LeDcUiou & Hen- 
riksen 2003). 
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merical convergence is reached when falls in this range of 
values. 

2) Any cold isotropic distribution of mass develops 
clumps at the onset of collapse. Roy & Perez (2004) seeded 
some of their initial mass distributions with homogeneous 
clumps (also spherically symmetric). Such initial conditions 
are difficult to duplicate for non-homogeneous spheres as 
done here owing to the tidal force of the background po- 
tential which does not define a spherical Roche boundary. 
For that reason, we limit our study to distributions of point 
sources with no internal degrees of freedom. This is only 
a minor setback however because clumps that form in our 
simulations during in-fall grow self-consistently from Pois- 
son (root-n) seeds. These small clumps merge as in-fall pro- 
ceeds which lead to the growth of a few large clumps just 
prior to maximum in-fall (see e.g. Fig. 7 of Aarseth, Lin & 
Papaloizou 1988). Thus the full process is highly anisotropic 
despite the choice of spherically symmetric initial conditions. 

3) Roy & Perez (2004) considered finite-Q (their param- 
eter 77) initial conditions only; we will show that the limit 
where Q — > gives rather different equilibria. In that sense, 
our survey complements theirs while focusing on more spe- 
cific initial conditions. 

It may be worth commenting that while we chose power- 
law initial conditions to generate power-law equilibria, an 
expectation drawn from past experiments with initial value 
problems of this kind, we found surprisingly that power-laws 
are recovered only for a sub-set of the parameter space, a 
result that turns on its head the long-held belief that the 
phase-space structure in equilibrium correlates strongly with 
the initial conditions chosen. While correlations are found, 
the end-results do differ significantly from those anticipated. 

All numerical studies using particle-based methods have 
used the virial ratio Q as a free parameter. This is partly 
justified on the grounds that there may not exist unique di- 
agnostics for the growth of ROI either in terms of a critical Q 
or a global anisotropy parameter in a spherical equilibrium 
(Palmer & Papaloizou 1987; Perez et al. 1996; see Merritt 
1999). Generally lower values of Q lead to deeper radial in- 
fall and more anisotropic velocity fields in equilibrium (more 
radial orbits), which favours ROfl Henriksen & Widrow 
(1997) have shown using a one-dimensional code that orbit- 
crossing which develops during in-fall leads to self-similar 
patterns which break off through a phase- mixing instability. 
Rapid phase-mixing softens the regime of violent relaxation 
that then sets in. Merral & Henriksen (2003) argue that this 
instability will develop more slowly when Q > 0. Conse- 
quently the similarity pattern persists longer in those cases. 
We will address this point in a study of three-dimensional 
'warm' collapses by comparing their outcome with those of 
cold collapses. 



^ The ROI is a Jeans type of instability developing in equilib- 
rium systems. The two-stream instability (in- and out-flow) may 
provide a more appropriate description of the outbreak of a bar 
during in-fall; see Barnes, Hut & Goodman (1986). 



2 Method 

2.1 Initial conditions 

2.1.1 Cold collapses 

Similarly to other authors (Aguilar & Merritt 1990 = 
A&M-I-90; Cannizzo & HoUister 1992 =C&H-|-92; Henriksen 
& Widrow 1997) we setup a series of spherically symmetric 
cold scale-free distributions of mass density 



po{r) oc {r/rsY 



(1) 



with ^5 7 < 5/2. The lower limit corresponds to homoge- 
neous spheres; the upper bound corresponds to systems with 
finite gravitational binding energy GM'^/r oc r'^-^^. Notice 
that all models with 7 > 1 have a diverging force field at the 
centre. A value of 7 = 3/2 was adopted as reference for the 
setup of the numerics and convergence of the parameters. 
The full range of models is listed in Table [1] 



2.1.2 Scaling the velocities with Q 

In order to assess how warm initial conditions infiuence equi- 
librium profiles, we ran also a few simulations in which the 
particles did not start from rest. Instead, they proceed from 
Dehnen (1993) density profiles and matching velocity field. 
The density profile of these models is given by 



p{r) 



(3 - 7) M 
47r 



ro 



(2) 



where 7 fixes the inner power-law index, (3 — 4 — ro and 
M are constants fixing the scales of length and total system 
mass. Because of known correlations between initial condi- 
tions and equilibria through incomplete relaxation (cf. §1.1), 
we expect initial- and equilibrium power-indices of the mass 
profiles near the centre to be equal. Of all the Dehnen mod- 
els the one with 7 = 3/2 gives the best match to a de Vau- 
couleur profile, as the reference power-index value adopted 
for cold collapses. 

We define the system virial ratio of kinetic to gravita- 
tional energy, Q, as 



Q 



2Ma^ 
W 



(3) 



where a is the three-dimensional velocity dispersion (we con- 
sider only irrotational models). Q = 1 for a virialised equi- 
librium, whereas Q < 1 ensures that all shells of constant 
mass converge to the barycentre of the system. Several stud- 
ies of gravitational collapse assign particle velocities ran- 
domly (e.g. McGlynn 1984; A&M-h90; C&H-F92; Boily et 
al. 1999). In the present study, however, we assign veloc- 
ities from the equilibrium distribution function as follows. 
The gravitational binding energy W is computed directly 
from ([2]) while the global mean square velocity dispersion 
is obtained from the Jeans equation in spherical symme- 
try (Binney & Tremaine 1987, §4.2). 

We then attribute velocities according to an isotropic 
Maxwellian velocity distribution constrained to satisfy lo- 
cally the first three moments of the Boltzmann equation. 
The velocities are then renormalised to achieve the desired 
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global virial ratio Q in (|3} . This approach has the advantage 
that the velocities are self-consistent with the mass profile of 
the system, and follows the strategy adopted by Barnes et al. 
(1986) in their study of the growth of ROFs in equilibrium 
systems. 

2.2 Choice of integration parameters 

2.2.1 Choice of units 

Power-law distributions were all truncated at a fixed radius 
rt = 2. We adopted units such that G = M = 1, which, 
together with rt, sets all scales in the problem. For an initial 
density profile p(r) = po{r/rs) the scales of length and 



of density po satisfy pov] - 
7)/47r. 

2.2.2 Time-step and free- 
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time 



The evolution of the A^-body models was followed on the 
Marseille Grape-5 systems (Kawai et al. 2000), using a spe- 
cially adapted TREE-code (Athanassoula et al. 1998). The 
code uses a fixed time-step which we set equal to 

St = tfi/3000 

where the nominal free-fall time, tg, is 



ts 



Stt 



32G(p) 



(4) 



where (p) — M / {A-Kr'^ / 'i) is the mean density. In computing 
units ts ~ 3.07(0). 

Direct application of the virial theorem (r ~ f't/'2) at 
constant M means that the system crossing time in equi- 
librium tcr = 2r/a\d = Ij-K ts — 0.637 ts will take a value 
close to half the free-fall time, or 1.95(6) A'^-body time units; 
we will use t^r = 2 A''-body time units for convenience when 
discussing the results. The code uses a Plummer smoothing 



GM 



to avoid divergences due to particle-particle interactions. 
The nature of the problem at hand requires a careful setup 
to ensure that global energy and angular momentum are 
preserved to good accuracy throughout evolution. The max- 
imum relative error 5E/ E (in percentage) measured during 
the complete evolution fluctuates between different runs, but 
remains of the order of, or better than, a few parts in a thou- 
sand for the adopted values of the opening angle, softening 
and time step. We discuss these values briefly. 

2.2.3 Choice of softening t 

As in all cases, the value of e should be tailored to the prob- 
lem at hand. The softening must be large enough to avoid 
collisions between particles, while being small enough to re- 
solve time-variations in the potential at all stages of evolu- 
tion. Our problem is particularly difficult, since it includes 
both very high density and very low density regions. Fur- 
thermore, particles do not stay all through the evolution in 



regions of similar density, but can visit regions of very dif- 
ferent density. 

We tried a number of softening lengths ranging from 
1/32 to 1/1024 and after many tests, some of which are 
described in section |4l we adopted a reference smoothing 
length of 1/512. The effect of changing this value and how 
this influences results are discussed in section IT3] 

Of course we can not trust our results for distances 
smaller than a few times the softening length. In practice 
we took a conservative radius of 5 x e as the inner-most 
radius to compute physical parameters (velocity dispersion, 
density, etc). 

2. 2. 4 Choice of opening angle 6^ 

The TREE method defines a critical angle 6c to control the 
accuracy of a limited expansion to the force field. The choice 
of 6c is a compromise between high accuracy and perfor- 
mance. In most simulations an opening angle 6c ~ 0.7 allows 
energy conservation typically to the order of 0.1% (Barnes 
& Hernquist 1996; Athanassoula et al. 2000). For collapse 
calculations, however, as for mergings, the time-variations 
of the potential are rapid and large, so such an opening an- 
gle may not necessarily suffice. We have thus carried out a 
series of calculations varying 6c and Plummer softening to 
determine which combination of parameters allows integra- 
tion to the desired accuracy. Our results proved robust for 
a wide range of opening angles (cf. Fig. [5]). We adopted a 
median value 6c ~ 0.4 for all our calculations. Simulations 
with an opening angle of 0.7 give a somewhat larger value 
of the energy variation. We thus did not include them in 
the analysis presented in the next sections. In general, runs 
that gave energy errors larger than 1% were discarded from 
analysis altogether. 



2.3 Physical description 

Starting with zero kinetic energy, all particles converge to- 
ward the centre of gravity which coincides with the centre 
of coordinates. Since the free-fall time 



ts oc 



(r/rs) 



7/2 



(5) 



increases with radius, several orbits cross at the centre be- 
fore a significant fraction of the mass has reached there. This 
feature allows us to treat the problem as an accretion prob- 
lem. At constant radius, the flux of matter through a shell 
during in-fall 



: 4-Kr^p{r,t)vr{t) cx r^-^^/^ oc t3(2-7)/7 



(6) 



increases with time when 7 2. The above relation was 
derived by following an in-falling shell of constant internal 
mass and so holds up to first orbit crossing. As orbits be- 
gin to cross at the centre, a pattern emerges from the origin 
and propagates outwards. At a given radius, out-going par- 
ticles eventually meet with in-falling material and the net 
mass in-flux rh drops. Thereafter a self-consistent equilib- 
rium is established over some « 20 A'^-body units of time 
(or, ~ 10 tcr); all our simulations ran for 80 A'^-body units 
of time, or ~ 40 dynamical times, to ensure stability of the 
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equilibrium configurations. Note that the classical two-body 
relaxation time at the system half-mass radius, trh, is given 
by (see Binney & Tremaine 1987) 

trh ^1 N 

1~ ~ 10 ln0.4iV 

so that for A'' = 100,000 particle calculations trh ~ 900 fcr 
is much larger than the total runtime of the simulations. In 
fact, the central relaxation time, defined in terms of the dy- 
namical time at the centre, can accommodate an tcr shorter 
by a factor 20, or an increase in density of 20^ = 400, and 
still remain clear of two-body relaxation effects. We checked 
explicitly for several cases that runs with this number of 
particles did not exceed a density contrast (central to half- 
mass values) of 100 (cf. Fig. [1} . We have performed a few 
runs with particle number A*' ~ 10*, to ease comparison 
with results from previous papers where this value of was 
used. Runs with A*" ~ 10* have trh ~ 160 tcr at the half 
mass radius; the maximum density contrast allowed is now 
~ 16, close to what was obtained (Fig. [TJ . These runs will 
likely have suffered a degree of two-body relaxation at the 
centre : to ensure that our analysis was not affected by re- 
laxation for this case, we also took snapshots at time t = AQ 
(half the evolution time) and verified that very similar con- 
clusions applied. Therefore only results obtained for t = 9>Q 
units of evolution will be presented. We note that the relax- 
ation time estimates are conservative since they do not take 
account of softening. This cuts off large-deflection angle en- 
counters and lengthens collisional relaxation (see e.g. Theis 
& Spurzem 1999; Athanassoula, Vozikis & Lambert 2001). 



3 Example : time-evolution of an 7 = 3/2 run 
3.1 Centre of density vs centre of mass 

Gravitational collapse leads to a spread of the individual 
particle energy range. Particles acquiring positive energy 
leave the system on a dynamical time-scale. Overall, on the 
order of 10% of the particles escape. The fraction of escapers 
is a function of both accretion index 7 and initial morphol- 
ogy. The combined effects may cause up to ~ 22% mass loss 
for homogeneous spherical systems (Table 1 of C&H-I-92; 
Boily 1994). As a result of the growth of radial orbit insta- 
bility and of anisotropic loss of unbound particles, the cen- 
tre of coordinates does not necessarily match the centre of 
mass of the bound particles. We found in most cases a small 
but non-zero linear momentum carried away by escapers. 
In all the analysis, and in particular when fitting isodensity 
curves and calculating the inertia tensors and their eigenval- 
ues, we recentered the coordinates to the center of density 
(i.e. the maximum density point) in order to avoid errors due 
to off-centered isophotes. The density maximum was identi- 
fied using the six-nearest neighbour scheme of Casertano & 
Hut (1985). The approach allows a much higher resolution of 
density profile around the centre, down to a few xe (Fig.[T|. 
The NEMO package utility functions greatly helped auto- 
mate the procedure (see |http : / /bima. astro.umd.edu /nemo / 1 
Teuben 1995). 



3.2 Inside-out morphology 

We investigated the morphology of equilibrium profiles in 
the centre of density frame in two different ways. First we 
sorted the particles in binding energy and removed the 20% 
particles with the largest energy (which includes all unbound 
particles). We then divided the remaining particles in four 
bins, each containing 20% of the initial particle number. In 
this manner we did not take into consideration unbound par- 
ticles and particles bound to the system but in a very-low 
density region: this imitates the effect of a tidal boundary. 
We also used an alternative approach, which consisted in 
sorting particles by increasing spherical radius and remov- 
ing the outermost 20%. In both cases we re-centered on the 
density maximum of the remaining particles. The two ap- 
proaches gave similar results. 

We computed the inertia tensor of selected particles, 
separately for each bin. The orientation of the system was 
set so that Cartesian axes matched the principal axes of 
the bound particles, with the positive x-axis coinciding with 
the semi-major axis. We then defined ellipsoidal rms axes 
a > b > c from the inertia tensor I as follows, for instance 
for the minor axis c 



where the (..) is the average of a quantity over A'^; particles 
of mass m in the i**" bin. From our simulations we recover 
spatially- and time-dependent axial ratios. 

Figure [2] graphs the time- evolution of the axis-ratios 
b/a and c/a for four bins for the duration of an A' = 800, 000 
run (d007 in Table [TJ. The configuration initially had a 
power index 7 = 3/2. Since we start from spherical sym- 
metry, both axial ratios are initially equal to 1. However, 
the rapid growth of radial motion leads to instability and 
the configuration quickly becomes triaxial (around t = 5 
on the figure). Axial evolution slows down rapidly, until 
after t ~ 20 when the curves flatten out and begin to fluc- 
tuate about average values. What morphological evolution 
remains is largely confined to the innermost 20% mass bin, 
seen both for the minor c-axis and the median b-axis. Thus 
the inner region continues to evolve, albeit slowly, within 
the overall relaxed structure. A close look at Fig. [2] shows a 
gentle but steady increase of the minor-axis ratio c/a, from 
« 0.38 to « 0.43, for the innermost three mass bins of the 
system. The minor-axis ratio of the outer-most particles 
show no indication of evolution. By contrast, the ratio b/a 
evolves significantly for the 20% most bound particles only, 
from 0.58 to ~ 0.65. This is suggestive of a trend toward 
axial symmetry as argued by Theis & Spurzem (1999) and 
Heller (1999) in their relaxation study of Plummer spheres 
(see also Curir & Diaferro 1994) . The example depicted on 
Fig. [2] suggests more rapid evolution in the inner, denser 
region of the system (shortest dynamical time-scale). We 
note that the initial conditions used here are without a 
harmonic (Plummer) core. The similar qualitative trends 
obtained from distinct sets of initial conditions supports 
the view that evolution toward axisymmetry is a generic 
feature of collision-less triaxial equilibria. 
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Figure 1. Density profiles of a series of runs with 7 = 3/2. The number of particles for eaeh case is indicated. The density was averaged 
over spherical shells containing 100 particles each, centered on the density maximum identified using the six-nearest neighbour scheme 
of Casertano &; Hut (1985). The larger simulations allow deeper probing of the inner region. Two runs with N = 100,000 particles but 
different random number seeds are displayed to illustrate scatter. The smoothing length and half-mass radius are marked with arrows, 
and a straight line of slope -3/2 (p oc r~^^'^) is shown for reference. 
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Figure 2. Time-evolution of the axis ratios for run d007. 



4 Morphology of the equihbrium profiles 

We now turn to a more quantitative inspection of our mod- 
els. A&M-I-90 found from their 5,000-particle runs that the 
bulk of their triaxial equilibria shaped by the radial orbit in- 
stability are prolate (see their Table 1). This result was con- 
firmed by C&H-I-92, who showed using 10,000-particle runs 
that the morphology of virialised objects varies in space, 
the equilibria being more triaxial prolate at the half-mass 
radius, and near-spherical in the inner and outer parts (see 
their Table 1 and Fig. 4). These studies set references against 
which to compare our results. In this section, we discuss the 
effect of varying numerical parameters (such as the number 
of particles and the softening) and model parameters (such 
as the index 7) on the equilibrium morphology; we consider 
parameterised fits to the density profiles in §5. 



4.1 Quantifying morphology 

A&M-I-90 introduced the parameter 



T = ^ (7) 
a — c 

to put the morphology of the relaxed equilibria on a quan- 
titative footing, r is bounded between (prolate spheroids) 
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Figure 3. Time-evolution of runs d007 (A^ = 800,000 particles, left-hand panel) and d022 (N = 200,000 particles, right-hand panel). 
We graph the morphological parameter r* defined in ([HJ for different 20%-mass shells, as indicated. 



and 1 (oblate spheroids). The degree of oblateness or pro- 
lateness varies within the bracket [0, 1] : configurations with 
r > i are oblate, otherwise they are prolate. 

In practice r is sensitive to round-off errors when a « 
6 ~ c. A clearer distinction between the three morphologies 
(triaxial, prolate and oblate) would be desirable. Consider 
instead the parameter 



which covers the range [—1, 1]; =0 for a spherical distri- 
bution, while r* > for 'oblate' triaxial distributions and 
< for 'prolate' triaxial distributions. The parameters r 
and r* are related to each other and both increase mono- 
tonically with oblateness, however note that r, allows to 
identify axisymmetric discs (r, = 1) or spindles (r, = — 1) 
unambiguously. For this reason we have chosen to use r, in 
our analysis. 

As an example, the evolution of r, with time is shown 
on Fig.O left-hand panel, for the 7 = 3/2 case displayed on 
Fig. [5] Although on the whole the system achieves oblate 
morphology rapidly, the aspect ratios and hence the diag- 
nostic T* remain strong functions of the volume sampled. 
The innermost 20% mass shell first becomes mildly prolate 
before shifting to oblate shaped (r« > 0) after some 30 time 
units of evolution. This clearly serves as warning against 
hastily classifying the system globally as either prolate or 
oblate. 

4.2 Dependence on A'^ 

We looked for trends in axis ratios c/a,b/a, and in r» as 
function of the number of particles used in the calculations. 
Recall that relatively large variations are expected for TV < 
10^ (BAK-l-02). To show that the results vary little whenever 
N > 100, 000, first we graph on the right panel of Fig.[3]the 



results of the same calculation as shown on the left but now 
with A'' = 200, 000 particles. The trends and values over time 
of r» are essentially the same as for the larger A*' — 800, 000 
calculation. This suggests that the numerics have converged 
for that range of values of A''. However, full convergence has 
to be demonstrated through comparison with a set of results 
for smaller particle number. 

We selected eight runs with A*' ranging from A*' = 10, 000 
to 800,000. Each run had a power-index 7 = 3/2. To min- 
imise root-A'^ fluctuations we averaged aspect ratios over 
time from 60 to 80 (end of the calculations) for a total of 
5 outputs. Table [5] lists their simulation parameters and as- 
pect ratios and r, achieved in equilibrium. (The aspect ra- 
tios were evaluated separately for four 20%-mass bins and 
for each output, however only values averaged over all mass 
bins are given.) We added to this compilation the results 
obtained by C&H-I-92 for this value of 7 (their index n) and 
N = 10,000, as well as results from A&M-I-90, who con- 
sidered cases with 7 = 1 and N = 5, 000. (A comparison 
of results from C&H-I-92 for the same 7=1 values with 
those of A&M-I-90 shows a good agreement between these 
two studies.) The values of r« computed from the data of 
A&M-I-90 were increased by 4-0.15 in account of the fact 
that they used a smaller 7 and in anticipation of the trend 
seen in plots of r» versus 7 (see §4.4 and Fig. [6] below). 

The results are graphed on Fig. |4l right-hand panel. 
We find predominantly oblate equilibria (r, > 0) for all our 
simulations with A'^ ^ 25, 000. Fig. 2] indicates less oblate 
structures from the A*' = 10, 000 cases than for larger- A'^ 
runs (open squares on the figure). We observe that the full 
range of axis ratios is reached only for N > 100, 000. Thus we 
find equilibria with ratios a /cm 2.6 for that order of particle 
number, well above the maximum a/c « 2 for 10,000-particle 
runs (cf. FigO middle panels). The same conclusion applies 
to the major-to-median axes ratio. 

The major- to minor-axis ratios of our 10,000 particle 
run, ranging from 1.6 to 2, all exceed those obtained by 
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C&H+92 for the same 7 = 3/2 case. What is more, the re- 
suhs we obtained for r* differ significantly from those of both 
A&M+90 and C&H+92 (we get r. « -0.308 from their Ta- 
ble f ): such results point systematically to prolate structures 
of equilibrium. Differences between our results and those of 
C&H-I-92 can not be attributed to particle number : below 
we make a detailed comparison of the initial conditions and 
software configurations to understand the origin of these dif- 
ferences. 



4.3 Smoothing length issues 

If e is large we expect a loss of resolution and a bias in the 
outcome of the calculations. Burkert (1990) and Boily et 
al. (1999) have discussed this in terms of a cut off in orbit 
deflection angle from the mean field and artificial saturation 
of the phase-space density. We looked for such biases in a 
series of calculations with varying smoothing length. 

Figures |4] and [5] (left-hand panels) graph the axis ra- 
tios and the averaged values of r, for the runs listed in 
Table O All runs had iV = 100, 000 particles. The aver- 
aging was done for the 80% most bound particles to facil- 
itate comparisons with C&H-I-92. The error bar displayed 
(of 0.22, bottom-left panel) is scatter. The trend seen in the 
data clearly indicates that higher resolution runs yield more 
oblate structures. Note the apparent convergence achieved 
for 1/e = 256, as r* and the axis ratios flatten out around 
that value. 

The results of Fig. |4] and [5] underline the importance 
of using both sufflcient number of particles and resolution 
for convergence of the equilibrium morphology. By contrast, 
the equilibria are not sensitive to the value chosen for the 
opening angle 6 (Fig. [S] right-most panels). If we compare 
values of r* at 1/e = 10 with those for 1/e = 512, we find 
a systematic increase of ~ 4-0.3, enough to bridge the gap 
between our results for A'^ = 10, 000 particles and those of 
C&H-I-92 seen on the right-hand panel of Fig. 4. Thus, the 
reduced resolution of their calculation accounts for the most 
part for the strong prolate morphology of their equilibria. 



4.4 Morphology in relation to 7 

The sensitivity of equilibrium parameters to numerical reso- 
lution suggested to us to seek a relation between equilib- 
rium morphology and initial density profile with higher- 
resolution calculations than was done in earlier work. On 
Fig. we graph the ratio of major- to shortest-axis a/c as 
function of the power index 7. Table |4] lists the parameters 
of the N = 100, 000 runs. The graph shows no trend with 
7, with values of a/c (here averaged over the 80% most- 
bound particles) ranging between 1.4 and 2. On the figure 
we added the linear fit to the data from C&H-I-92 for the 
same quantity for comparison. The results of an A = 10, 000 
and an A = 25, 000 calculation, with the same linear res- 
olution e = 1/512, are also added in Fig. [la,. The 10,000- 
particle result falls close to the straight line from C&II-I-92, 
while the 25,000-particle run gives a ratio comparable to 
those obtained with 100,000-particle calculations. This con- 
firms our expectation that 10,000-particle calculations do 
not quite reach maximum potential during violent relaxation 
(cf. BAK4-02) and as a result do not achieve in equilibrium 



the same final (converged) morphology of larger-A calcula- 
tions. Note that the axis ratios for the cases 7 = and 7 = 2 
differ significantly from those obtained with other power in- 
dices. Nevertheless, the full range displayed on the y-axis of 
Fig. |6^ brackets E2-E5 morphological types, which are by 
no means exceptional values when compared with those of 
elliptical galaxies. 

On Fig.lBJj we plot the values of r, as a function of initial 
power index 7. As before, the systems were split in equal 
20%-mass shells to compute r, at times ranging from 60 to 
80 A-body units. The scatter seen on Fig. results form 
fiuctuations between different mass bins. The solid line and 
filled symbols indicate averages. In contrast to the results 
for the aspect ratios, we find a well-defined trend of the 
parameter r* versus power index 7 (Fig. ^p) such that r* 
increases on the mean with increasing 7, and r* > for 
7 ^ 5/4. The case with 7 = 2 yields the most oblate and 
axisymmetric equilibrium of all the cases we have explored. 
Data points for 7 = 1 and lower are only mildly oblate on 
the mean, while the case 7 = 1/2 in fact gave a mildly 
prolate structure. However, the morphology of all runs with 
7 < 3/4 is consistent with axisymmetry or perfect-ellipsoid 
morphology (r* = 0), owing to the scatter. Seen in this light, 
we would conclude that the morphology of systems with 
7 ^ 3/2 is oblate for all mass shells, whereas those with 
lower 7 values are a mixture of oblate and prolate shells. 



5 Global and inner density profiles 

Fig. [1] displays typical spherically-averaged density profiles 
for our simulations with 7 = 3/2. Our basic expectation 
from incomplete relaxation is that the density profile will 
match the power-law of the initial conditions on small scales. 
Power-law asymptotes at small and large radii are a generic 
feature of dissipation-less relaxation in cosmology (Navarro 
et al. 1997; Moore et al. 1998). It is natural to expect them 
also in the present, more restraint, context. We therefore set 
out to fit the radial density profiles of our equilibria with 
a continuous function of radius. All density profiles where 
averaged over spherical shells containing 100 particles. We 
varied the number of particles per shell (200, 500) but this 
did not lead to significant improvements, save for the very- 
large A d007 model where 500 particles still allowed to probe 
a sufficient range in radius while keeping the noise level low. 

C&H-I-92 had found a good fit to their equilibrium 7 = 1 
run averaged over spherical shells using an Hernquist (1990) 
density distribution, 

p{-r) = 7, 7 — — (9) 

ZTT r (ro + rj'^ 

where M is the total system mass, and ro a free length. We 
also fitted the equilibrium profile obtained from an 100,000- 
particle 7 = 1 run (d013, cf. Tabled with Q. We first 
evaluated the density locally following the scheme of Caser- 
tano & Hut (1985) as coded in the NEMO analysis package. 
The rms axis ratios were found as in §3, again covering the 
innermost 80,000 particles only. This allowed us to round up 
the distribution, by shifting the particles in x-y-z inversely 
as the axial ratios obtained : thus the triaxial configuration 
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Figure 4. Global morphology parameter t* as function of inverse smoothing length 1/e (left-hand panel) and particle number N (right- 
hand panel). The vertical tick mark on the left panel is an error estimate derived from the scatter seen in the axial ratios; the data 
points are mean values for the innermost 80% mass. The data for the 5,000-particle runs are taken from A&;M-|-90 but shifted upwards 
by -1-0.15 (see text for details). The symbol 'f is the mean value lifted from C&H+92 for N = 10, 000. 
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Figure 5. Axis ratios a/c and b/c as function of the opening angle 9 (right-hand panels), particle number N (middle panels) and 
resolution 1/e (left-hand panels). Simulation parameters are listed in Tables [2l and l3l Values obtained for each of four 20% mass bins are 
displayed with different symbols. 
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Figure 6. (A, left-hand panel) Ratio of major axis a to shortest axis c for runs with A'^ = 100,000 particles (a) and different power 
index 7. Results with 10,000 (□) and 25,000 (O) particles are also displayed for the case 7 = 3/2. The straight dotted line is the fit 
a/c = 7/5 + 1.3 lifted from C&H+92. (B, right-hand panel) Shape parameter r* versus initial power index 7. The horizontal line indicates 
the boundary between oblate (positive) and prolate shapes. On the mean the structures are mostly oblate, more so for 7 > 1. Averaged 
values for each 7 are linked with solid lines. The star marks the average result of 5,000 particle 7 = 1 runs from A&:M-|-90; the dagger 
at 7 = 3/2 marks the average result of 10,000-particle runs by C&H-I-92. 




Figure 7. (A, left-hand panel) The equilibrium density profile of run d013 versus semi-major axis length. The dashed line is a fit of the 
data with an Hcrnquist model (Eq. [9]l. The straight line in the upper-left has a slope of -1 and the half mass radius is given by a 
vertical tick mark. (B, right-hand panel) The logarithmic derivative of the density profile shown on panel (A). The solid line is obtained 
from Eq.[9] The vertical tick marks indicate (from left to right) the 1%, 10% and 20% most bound mass fraction. 



was straightforwardly transformed to a sphere of the same 
volume. 

The result is shown on Fig. [7^. The curve shown has 
ro = 0.288 with M — 1. Note that the half-mass radius 
(rh — 0.53) obtained numerically is significantly less than 
the Hernquist function = (1 + \/2) '"0 = 0.69(5) owing to 



the mass deficit at large distance. Inside r ~ 10 the mass 
integrated from ([9| is 92.5% of the total, while the mass 
found numerically reaches 91.2% of the total. However , the 
data and analytic curve are in complete disagreement for 
r > 10. Nevertheless a comparison of this fit to the one 
shown on Fig. 8 of C&H-I-92 indicates that our larger-A?^ 
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Figure 8. The logarithmic derivatives for equilibria obtained from different initial power-law index ■y. N = 100, 000 in all cases but for 
the case 7 = 3/2 which had N = 800,000 particles. The vertical tick marks indicate (from left to right) the 1%, 10% and 20% most 
bound mass fraction. 
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Figure 8. - continued 



simulations give an ever closer agreement to the Hernquist 
profile. The fact that both profiles are relatively well fitted 
by the Hernquist law is unexpected, since the morphology of 
the two equilibria are very different from each other (oblate 
here, cf. Fig. [BJj, and prolate in C&H+92). 



5.1 Logarithmic derivative : 7 = 1 

One reason why a global function does a poor job of fitting 
the mass distribution is the irregular, noisy profiling outside 
(about) the half-mass radius. We illustrate this with a series 
of graphs of the logarithmic derivative of the density as func- 
tion of radius. A piecewise linear fit to the finite-difference 
scheme was done, such that 



dlogr . log(ri/ri+fc) 



(10) 



where i := [1, 1000 — k] is the index of a mass shell, and 
fc « 30 is constant. The value of k is chosen so as to sample a 
sufficiently large radial increment and avoid round-off errors 
in (|10p . We found that k could be increased to ~ 200 without 
affecting the overall picture much. We then used a least- 
square fitting routine (Press et al. 1992) to find the best 
linear fit over k points. The result was centered on the mean 
radius in the interval + k. 

On Fig. [Tja we graph this derivative for the run 7 = 1 
listed in Table [l] We have indicated the radii of the inner 1%, 
10% and 20% mass shells with vertical ticks on the figure. It 
is clear that the derivative is smooth and well fitted up to the 
20% mass shell; however beyond that point large fluctuations 
are seen which indicate grainy or irregular substructures. 
Note also that the data match very well the logarithmic 
derivative of (|9]), given by a solid line on the figure, from 
the innermost 0.1% mass bin up to roughly the 20% mass 
radius. Note the apparent clustering around dlogp/dlogr ^ 
—2.2 also between the 20% and half-mass shells (rj, ~ 0.53). 



Beyond the half-mass radius the derivative oscillates wildly 
between values of —3 and —4; beyond r = 1 or so the data 
is more erratic. 



5.2 Logarithmic derivative for all 7's 

It may come as a disappointment that the logarithmic 
derivative on Fig. [Tja has not converged to the anticipated 
constant value = — 1, despite the very small innermost mass 
fraction (0.1%) sampled. We asked whether the same was 
true of all the runs with different 7's. We therefore repeated 
the procedure for all runs listed in Table |4l except for the 
case of 7 = 3/2 where we took the results from the 800,000 
particle run. The results are displayed on Fig. [S] Here again 
the vertical ticks indicate 1%, 10% and 20% mass shells. For 
small values of 7, dlogp/dlogr does not converge to a fiat 
value as we approach the center. However, as we consider 
larger values of 7 we find hints of such a convergence, par- 
ticularly for 7 > 3/2. Once more we find a suggestion of a 
leveling off around dlog p/dlog r ~ —2.2 at radii close to the 
20% mass shell for all the cases with 7 ^ 1. In their treat- 
ment of the Boltzmann equation in one dimension, Hozumi 
et al. (2000) also found a power-law fit of index ~ —2.1 
around the half-mass radius of their systems. 



5.3 Circular velocity 

On Fig. |9] we graph the (less noisy) circular velocity 
Vc = \/ GM/r as function of the radius and of the inte- 
grated mass. Note that Vc — > constant when p(r) oc r~ . On 
the figure we also added a straight line (oc r or oc 
that Vc would follow if the density flattened out at the 
centre. For our simulations, Vc rises slowly as a function of 
radius and we find Vc ~ r''''^^ gives a rough fit for runs with 
7 < 1. Starting with the curve 7 = 1, one can see a plateau 
appearing at a mass fraction M{< r)/M « 0.1, which 
becomes more prominent for increasing 7. The circular 
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Figure 9. Rotational velocity Vc plotted for different 7 as function of radius (left-hand panels) and of integrated mass M(< r) (right- 
hand panels). Lower set : < 7 ^ 1. Higher set : 1 < 7 ^ 2. The straight lines are oc r (left-hand panels) or oc M^^^ (right-hand panels) 
corresponding to a constant mass density. 



velocity drops off in a Keplerian tail at large radii, in all the 
cases. We remark that the maximum of Vc nearly always 
occurs at the mass fraction M{< r)/M ~ 1/3 while the 
maximum itself is a non-monotonic but nearly flat function 
of 7 (save 7 = 0, see Fig. llOp . Therefore, while none of the 
equilibria can be mapped into another (different rotation 
curves, mass proflles), the results of collapse calculations 
for different 7's (and hence the accretion rates or history) 
are fairly similar in terms of max(i)c), and of corresponding 
mass fraction and radius, a situation that parallels current 
debates about the universality of resolved circular velocities 
and virial masses in cosmology (see e.g. Power et al. 2003; 
Navarro et al. 2004). Only the two cases of 7 = and 2 
stand out. These cases are particularly interesting, since 
they produce the largest and second largest maximum Vc, 
but the smallest corresponding radius and mass fraction 
(bottom panel. Fig. [111)1 . 



5.4 Power-law limit 

The search for a flat logarithmic derivative can be done 
through a functional fit which admits a linear regime in log r. 
We modified ([SJ slightly to look for fits to the derivatives of 
the form 

_ dlogp r'^ 

j{r) = — = -7/ -72logr-/3— - — - (11) 

d log r + 

where ro and (3 are as before, the symbol 7/ is used to dis- 
tinguish the final and initial values of 7 and the constant 
72 is equal to the second derivative at small radii to O(r^). 
We imposed /3 = 4 — 7/ so three parameters remain free. A 
strict power-law at short distances must yield 72 = 0. We 
applied the functional j{r) to three cases with 7 = 3/2, 7/4 
and 2. The results are shown on Fig. 1111 The data points 
are those of Fig. [8] in each case, to which we have added 
j{r) (solid line) and j'{r) = dj{r)/dlogr (dashed line). We 
find j'{r) < everywhere for 7 = 3/2, however a fit with 
72 = gave sensible results for both 7 — 7/4 and 2. Even 
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Figure 10. Maximum of tlic circular velocity Vc as function of 
the initial power- index 7 (top panel). The lower panel shows the 
mass fraction (m) and radius (rm) at maximum. 



SO the power- law regime in the case 7 = 7/4 does not cover 
more than a few percent of the system mass, while it extends 
roughly out to 10% of the mass for 7 = 2. 

Although the anticipated relation 7/ = 
lirnr^o —dlog p/dlogr ^ j = constant at small radii 
is not recovered, a strong correlation is found between 7 
and the minimum of — dlog p/dlog r read off Fig. [S] We 
graph — cJlogp/dlogr for the innermost 1% (solid) and 
0.1% (dotted) mass shells on Fig. 1131 The diagonal on that 
figure is the sought equality 7/ = 7. We find a clear trend, 
especially for 7 ^ 3/2, the strongest for the innermost 0.1% 
of the mass. For runs with 7 < 1 the correlation is more 
suggestive, even for the innermost 0.1% mass shell. 

As long as we are concerned with the radial profile of 
the systems at small r, these results show that memory of 
the initial configuration, while still there, involves only a 
very small fraction of the total mass. Correlations in binding 
energy have been noted by several authors (see e.g. van Al- 
bada 1982; Henriksen & Widrow 1999), so the initially most- 
bound particles are those contributing to correlations seen 
on Fig. 1131 We deduce that well-resolved equilibria could not 
have inner cusps with a logarithmic slope lower than 7 of 
the initial conditions. But significant deviations already at 
the 1% mass fraction means that this conclusion effectively 
applies to a very small volume. We may ask whether the 
velocity field also shows traces of the collapse, and in what 
measure. 



6 Phase-mixing and relaxation 
6.1 Cold initial conditions 

Having brought to light a trend for the run of density in equi- 
librium as function of the initial power index 7 (Fig. I13|l , we 



now ask whether a similar relation exists for the velocity 
field. Repeated exchanges of kinetic energy modify a parti- 
cle's velocity by 5v (say) each time. When these boosts in 
velocity are randomly oriented the ensuing velocity vector 
is uncorrelated with the initial v. This will be true in an 
average sense when the sum of all perturbations exceeds the 
norm of the velocity vector, i.e. ^^Sv^ > v'^ . For systems 
collapsing from rest, v = 0, this condition will be met first 
by particles that do not acquire large velocities while the 
background potential changes rapidly. From these consider- 
ations, we would expect the velocity field to be the one aris- 
ing from violent relaxation for all orbits with relatively small 
in-fall velocities. Because of the enforced spherical symme- 
try of the initial conditions, all orbits are radial {v = Vrv) 
at the on-set of collapse. The work by the gravitational force 
Svr oc rWr4> '"'^"^ ^ for all 7's < 2. The kinetic energy 
remains small for all particles orbiting near the centre and 
hence small fiuctuations in energy will wipe out memory of 
the coherent radial collapse. 

A Maxwellian coarse-grained velocity distribution func- 
tion (d.f.) (more precisely: a sum of Maxwellians) is the sig- 
nature of violent relaxation, a fully stochastic redistribution 
of kinetic energy. Merral & Henriksen (2003) and Iguchi et 
al. (2005) have obtained good agreement with this predic- 
tion by integrating the collisions-less Boltzmann equation 
directly, and evaluating the velocity d.f. f{y) in equilibrium 
at the origin of the coordinates. The same procedure can 
not be done here due to finite resolution. Instead, we eval- 
uated f{v) for an 100,000-particle 7 = 3/2 calculation by 
sorting and binning in velocity space all particles within a 
given mass fraction. We picked a fraction M{< r) /M = 2.5% 
which samples the regime where the logarithmic derivative 
is (roughly) linear with the logarithm of the radius, that is, 
the density p oc r*^'""" with k some numerical constant (cf. 
Fig[8]) and much larger than the softening length e ~ 0.002. 
For comparison, we also computed f{y) for all particles, tak- 
ing care to normalise the d.f. so that it integrates to unity 
in each case. 

The results are displayed on Fig |12l We find excellent 
agreement with a Maxwellian profile for the inner 2.5% 
mass saniple despite the rough dataset used (2,500 parti- 
cles onlyjj. Applied to the system as a whole, we find an 
overabundance of low- velocity particles, while overall the 
fit is not as good (see Fig. 1121 right-hand panel) . The ex- 
cess probability density at low- velocities with respect to an 
isotropic Maxwellian velocity d.f. may be understood as an 
overabundance of stars on low-energy radial orbits. Stars on 
such orbits spend more time at apogalacticon, where their 
velocity is relatively low. This is confirmed by looking at 
runs of the isotropy parameter /?* = 1 — 1/2 {a^/ a^) dis- 
tinguishing between perpendicular (_L) and radial velocity 
components (top panels. Fig. I12p . Clearly /3, = for an 
isotropic velocity field, and /3« > when the system shows 
an excess of radial motion. We find the inner region more 
nearly isotropic (left-hand set), while globally /3» assumes 
larger values and increases as decreases (right-hand set). 

In conclusion, we find a good fit to the inner region 



^ A repeat with the 800,000-particle run, or by averaging over 
time intervals 3> the local dynamical time, only confirms the qual- 
ity of the agreement. 
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Figure 11. Logarithmic first and second derivatives versus radius for the three cases with the largest initial power index 7. The solid 
line is a fit (j) and the dashed line the derivative of that fit {j', top-most curve). A strict power-law is ruled out in the case 7 = 3/2 
everywhere. For 7 = 7/4 a power-law inner profile out to fa a few percent of the mass is representative; while for 7 = 2 the constant 
power-law region reaches fa 10% of the total mass. 



velocity d.f. with a Maxwellian, while globally this does not 
hold. These results go in the same directions as the findings 
by others (e.g., Funato et al. 1992; Merrall & Henriksen 
2003; see also Bertin & Trenti 2003) but for fully three- 
dimensional calculations. By the time particles fall in with 
large velocities, much orbit crossing (and therefore phase 
mixing) has gone on in the centre. As a result, the potential 
varies more smoothly in time and these particles preserve 
significant radial motion. This is seen in animated snapshots, 
when the last mass shell has collapsed and particles rebound 
from the central region. A fraction of them are unbound and 
leave the system. 

6.2 Warm initial conditions 

We wish to assess how equilibrium profiles of simulations 
with Q 7^ initially differ from those starting from rest. 
Since kinetic energy provides pressure support, a run with 
Q > will collapse more slowly (reduced mass in-fall, rh{t)). 
The expectation drawn from simulations with different 7's is 
that a reduced rate of in-fall (large 7, cf. Eq.|6]) would favour 
oblate morphology in equilibrium (cf. Fig. [6]3). A&M-I-90 
have shown that sufficiently warm collapse simulations shut 
off ROI's and preserve the symmetry of the initial config- 
uration. (The limiting case of an initial Q — 1 equilibrium 
must trivially preserve spherical symmetry.) The question at 
stake here is whether ROI modes of instability are washed 
out progressively as Q increases from zero, leading to more 
and more spherically symmetric equilibria; or whether more 
oblate equilibria first develop for Q > but small, before 
being washed out. 

In order to answer this question we ran a series of sim- 
ulations with warm initial conditions. We set up Dehnen 
(1993) density profiles, Eq. with 7 = 3/2, /3 = 5/2, 
ro = 0.756.., which are truncated at radius rt ~ 3.76 so that 
the density profiles of these models match the density of 
the power-law models at the centre. This choice of scales al- 
lowed us to use the same numerical setup (smoothing length 
and time-steps) as for the power-law computations since the 
dynamics at the centre is identical in both cases. The to- 
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Figure 13. Value of the logarithmic derivative —ci log p/d log r as 
function of the initial conditions power index 7 for two mass shells 
of 1% (solid line) and 0.1% (dash). The latter is the minimum 
positive value obtained for each case. The diagonal broken line is 
equality. 



tal mass M(< rt) = 1 as before. As explained in §2, we 
attributed velocities according to an isotropic Maxwellian 
velocity d.f. constrained to satisfy locally the first moments 
of the Jeans equations (Binney & Tremaine 1987) . The out- 
come of such studies are not sensitive to details of the ve- 
locity field but to the global value of Q (Barnes et al. 1986; 
A&M-f90). 

6. 2. 1 Similarity patterns 

A self-gravitating collapse proceeding from power-law ini- 
tial conditions soon develops single-valued orbital patterns, 
such that any star follows the same orbit as any other 
once the scales of mass and lengths are renormalised un- 
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Figure 12. The velocity d.f. f{v) as function of v'^ . Note tlie logarithmic scale. The function was constructed for the innermost 2.5% 
particles at the end of the simulation (left-hand panel) . The same quantity computed for the system as a whole is shown for comparison 
(right-hand panel). The dashed curve is the same Maxwellian on both panels; the axes were rescaled for better display. The top panels 
show runs of the anisotropy parameter /3* = 1 — 0.5 (cr^/tr^); /3« > for systems with excess radial motion. The inner part (left-hand set) 
is more isotropic, while globally the system shows an increasingly anisotropic velocity field as we reach the edge of the system {v^ — ► 0). 



der a time-transformation. The trajectories of the stars map 
out a unique path and the flow is said to be self-similar. 
The growth of similarity patterns saturates to give way to 
a phase-mixing instability leading to virialised equilibrium 
(Henriksen & Widrow 1997, 1999). An example of this pro- 
cess is displayed on Fig. 1141 which compares the run of the 
radial velocity with radius for two time sequences, one corre- 
sponding to a cold collapse (d002) and the other to a warm 
one (h002). From this, and other similar plots for other sim- 
ulations and other times (not shown here) , it is clear that the 
self-similar pattern remains well defined for a longer time in 
the warm run than in the cold one. This is in good agreement 
with Merral & Henriksen (2003) , who argued that warm ini- 
tial conditions slow down the transition to a phase-mixing 
instability. This trend, however, is not monotonic. Indeed, 
for very hot runs, such similarity patterns hardly show up. 
Thus, intermediate values of Q, of the order of 0.1 or 0.2, 
are the optimum for their formation and maintenance. 



We also note that, for the warm simulation, the radial 
velocity dispersion is relatively large at the earliest time, and 
for sufHciently high Q values some particles even have Hr > 
and outward motion. As time runs, the radial flow becomes 
cooler (the one-stream at large radii becomes thinner). In 
the late stages the one-stream in-falling material is more 
and more fine-tuned (single- valued), as in the cold collapse. 
This observation will serve us below when interpreting the 
structural evolution of the models during collapse. 



6.2.2 Oblate, and then not 

Table [5] lists the overall morphology in equilibrium of a 
series of Dehnen models with 7 = 3/2 and different Q 
initially. The total runtime was 80 time units in each 
case. All collapses with Q > 0.4 (not listed in Table [Sjl 
led to highly spherical equilibria. For Q ^ 0.4 we get 
either prolate or oblate morphologies. For the coldest cases 
{Q ^ 0.05) the equilibria are oblate, as the scale-free mod- 
els. With increasing Q, the equilibria first shift to prolate 
and then to spherical symmetry. This is rather different 
from the monotonic trend we would have anticipated, 
had the aspherical modes of instability been washed out 
progressively with increasing Q. The transition is very 
sharp around Q ~ 1/10. Two possible explanations of this 
phenomenon are on hand. The first invokes the growth of 
a bending-mode of instability, similarly to self-gravitating 
thin discs (see Merritt 1999). Here, the highly elongated 
prolate morphology of the cold (low Q's) runs implies 
motion mainly down the semi-major axis. As Q is reduced 
and the bar achieves a highly prolate shape, an off-axis 
mode of instability will develop owing to the large velocities 
of the particles parallel to the semi-major axis. Comparing 
centrifugal and restoring forces in the case of a thin bar 
distorted by a sinusoidal bend of wave- number k = 27r/A 
and amplitude A immediately gives the condition Aft > 
constant for growth of the mode (the precise value of the 
constant is of no interest here). A similar relation applies for 
discs. Clumping modes of short wavelengths (large k) will 
develop more easily in cooler systems, whence the observed 
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Figure 14. Phase-space sections showing the radial velocity as function of radius during in-fall for two runs, one starting from zero 
velocities (d002, left-hand set) and the other from sub-virial velocities {h002, right-hand set). A similarity pattern shows up in the form 
of concentric rings, which are progressively deleted at later times. 



transition to oblate morphology. An alternative explanation 
is that a two-stream mode of instability develops more fully 
in the deep potential well of the colder collapses, when 
the in-falling material reaches larger velocities and more 
orbit-crossing takes place at the centre. Clearly a stability 
analysis well into the non-linear regime is required to de- 
termine which one of the two (or other) types of instability 
prevails. Such an analysis goes far beyond the objectives 
of the present paper and will not be attempted here. The 
close link between equilibrium properties and initial virial 
ratio was noted long ago for spherically averaged values 
(van Albada 1982; Mclynn 1984; A&M-I-90), however, to 
the best of our knowledge, the transition from oblate to 
prolate morphology and then to spherical symmetry as Q 
increases has not been stressed before. 

We close this section with a comparison of the cold 
Q — collapse from a power-law profile with the corre- 
sponding Dehnen model. The cold Dehnen run listed in Ta- 



ble [S] reached r* « 0.4 in equilibrium, i.e. more oblate than 
what was obtained from an 7 = 3/2 power- law initial con- 
ditions (cf. Fig. |4] and Table |4ll . Since the mass profile of 
a Dehnen model is steeper at large radii than p oc r~^^'^, 
mass shells falling in from large r take longer to reach the 
origin in that case. In other words, the rate of mass in-fall of 
the Dehnen sphere is lower compared to the case of collapse 
from a power-law profile. This further supports our claim 
that a reduced rate of mass in-fall rh (cf. Eq. |6]) leads to 
more oblate equilibria, at least when the initial virial ratio 
Q = Q. 



7 Summary and discussion 

Using AT-body computer simulations, we have shown that 
collapsing self-gravitating spheres develop oblate or prolate 
triaxial figures of equilibrium. In a study of numerical con- 
vergence, we have shown that an insufficient number of par- 
ticles N , or linear resolution e, can give wrong results for the 
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morphology. As A*' increases we find evidence for a slow drift 
from prolate to oblate morphology; the trend is similar but 
more pronounced when we reduce the smoothing length e. 
Previous studies had found prominently prolate structures 
of equilibrium, but this appears to have been due partly to 
the numerical setup used. The degree of symmetry of the 
virialised structures is highly sensitive to the phase-mixing 
which takes place during in-fall. We used scale-free power- 
law initial density profiles and found that steeper powers 
lead to more oblate equilibria. We also noted that increasing 
the initial virial ratio Q > leads to prolate equilibria. In- 
creasing Q to yet larger values has the effect of maintaining 
spherical symmetry, by shutting down aspherical modes of 
instability (Barnes et al. 1986). For very low-Q calculations, 
we noted that the prolate equilibria give way to oblate shape 
and have invoked two mechanisms that likely play a role in 
this transition (§6.2). A full analysis of that phenomenon is 
deferred to a future study. 

All our simulations lead to peaked central regions. 
When fitting the run of spherically-averaged volume den- 
sity with radius we found that the logarithmic derivative of 
the profiles converges only slowly to a power-law, at least 
when 7 3/2 (Fig. |8}. The logarithmic derivative in the 
central region of the equilibrium systems is well correlated 
with 7 when 7 3/2 (Fig. 1131) : in all cases, the central cusps 
would include only a very small fraction of the total mass, 
on a scale where the impact of baryonic (i.e., gas) physics 
would not be negligible. CoUisional relaxation effects are al- 
ways a worry in dense central regions. We checked that the 
collisional time of our simulations is too long, even at the 
centre, for two-body encounters to play a role. Power et al. 
(2003) conducted a thorough analysis of these eflects. We 
find our simulations on the safe side of their resolution cri- 
terion (see their Fig. 14). This boosts our confidence that 
the properties of the cusps and the evolution we observed at 
the innermost 20% mass shell is a genuine effect of the com- 
plex triaxial structure of the equilibria, and not a numerical 
artifact due to two-body relaxation. 

We also found that steep initial density profiles lead 
to equilibria with a leveling off in the circular velocity 
as function of radius. This suggests a run of density in 
equilibrium p oc for a non- negligible fraction of the 
mass. This feature was also reported by Hozumi et al. 
(2000), who approached the problem from the angle of the 
collision- less Boltzmann equation. 

We have observed evolution in time for the innermost 
region of the equilibrium profiles, such that the small-scale 
profile shifts away from triaxiality, and toward oblate axi- 
symmetry (see also Theis & Spurzem 1999; Heller 1999). 
This would certainly impact on the orbital structure of the 
equilibrium. In a follow-up investigation, we will explore the 
orbital structure of the equilibria obtained here in details 
to bridge over with models of triaxial ellipticals based on 
distribution functions. 
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Table 1. List of runs with cold scale-free initial conditions. 



Name 


N 




lA 


7 




[103] 






dOOl 


100 


0.7 


512 


3/2 


d004 


100 


0.6 


512 


3/2 


d003 


100 


0.5 


512 


3/2 


d002 


100 


0.4 


512 


3/2 


d005 


100 


0.3 


512 


3/2 


d006 


100 


0.2 


512 


3/2 


d019 


10 


0.4 


512 


3/2 


dOlO 


25 


0.4 


512 


3/2 


d008 


32 


0.4 


512 


3/2 


d009 


50 


0.4 


512 


3/2 


d020 


100 


0.4 


512 


3/2 


d022 


200 


0.4 


512 


3/2 


d023 


400 


0.4 


512 


3/2 


d007 


800 


0.4 


512 


3/2 


dOll 


100 


0.4 


512 





d012 


100 


0.4 


512 


1/2 


d015 


100 


0.4 


512 


3/4 


d013 


100 


0.4 


512 


1 


d016 


100 


0.4 


512 


5/4 


d017 


100 


0.4 


512 


7/4 


d014 


100 


0.4 


512 


2 


d031 


100 


0.4 


32 


3/2 


d030 


100 


0.4 


64 


3/2 


d025 


100 


0.4 


128 


3/2 


d024 


100 


0.4 


256 


3/2 


d021 


100 


0.4 


512 


3/2 


d026 


100 


0.4 


1024 


3/2 



Table 2. Morphology of equilibria obtained for runs with different particle number A'^. The smoothing is e = 1/512 and 7 = 3/2 in all 
cases. Values are averages over the innermost 80% system mass. 



Name 


N 


h/a 


c/a 


< T, > 




[10^] 








d007 


800 


0.76 


0.45 


0.21 


d023 


400 


0.81 


0.47 


0.28 


d022 


200 


0.76 


0.47 


0.20 


d002 


100 


0.78 


0.45 


0.27 


d009 


50 


0.72 


0.48 


0.06 


d008 


32 


0.77 


0.52 


0.12 


dOlO 


25 


0.78 


0.52 


0.12 


d019 


10 


0.80 


0.57 


0.10 



Table 3. Morphology of equilibria obtained for runs with different smoothing e. The number of particles is = 100,000 and 7 = 3/2 
in all cases. Values are averages over the innermost 80% system mass. 



Name 


1/^ 


b/a 


c/a 


< r* > 


d031 


32 


0.65 


0.43 


-0.03 


d030 


64 


0.70 


0.44 


0.08 


d025 


128 


0.72 


0.44 


0.14 


d024 


256 


0.86 


0.47 


0.38 


d021 


512 


0.70 


0.40 


0.17 


d002 


512 


0.82 


0.55 


0.30 


d026 


1024 


0.98 


0.56 


0.47 



On the equilibrium morphology of systems drawn from spherical collapse experiments 21 



Table 4. Morphology of equilibria obtained for runs with different initial power index 7. The smoothing is e = 1/512 and A'^ = 100, 000 
in each case. Values are averages over the innermost 80% system mass. 



Name 


7 


b/a 


c/a 


<T, > 


d012 


1/2 


0.65 


0.48 


-0.13 


d015 


3/4 


0.64 


0.44 


+0.08 


d013 


1 


0.74 


0.45 


+0.14 


d016 


5/4 


0.71 


0.44 


+0.10 


d002 


3/2 


0.78 


0.45 


+0.27 


d017 


7/4 


0.88 


0.49 


+0.39 


d014 


2 


0.98 


0.58 


+0.47 



Table 5. Parameters of the calculations based on 7 = 3/2 Dehnen models (Eq. 2). The smoothing length is e = 1/512 and = 100,000 
in each case. 



Name 


Q 


b/a 


c/a 


< T, > 


Comment 


hOOl 


0.0 


0.98 


0.61 


0.42 


Highly oblate 


h009 


0.05 


0.63 


0.52 


-0.18 




h002 


0.1 


0.63 


0.62 


-0.41 


Highly prolate 


h008 


0.15 


0.74 


0.73 


-0.21 




h003 


0.2 


0.98 


0.96 


0.00 


(quasi-) spherical 


h004 


0.4 


> 0.98 


> 0.98 


0.00 


spherical 



